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Abstract 

Deviations of grid frequency from the nominal frequency 
are an indicator of the global imbalance between genera- 
tion and load. Two types of control, a distributed propor- 
tional control and a centralized integral control, are cur- 
rently used to keep frequency deviations small. Although 
generation-load imbalance can be very localized, both 
controls primarily rely on frequency deviation as their in- 
put. The time scales of control require the outputs of the 
centralized integral control to be communicated to distant 
generators every few seconds. We reconsider this con- 
trol/communication architecture and suggest a hybrid ap- 
proach that utilizes parameterized feedback policies that 
can be implemented in a fully distributed manner because 
the inputs to these policies are local observables at each 
generator. Using an ensemble of forecasts of load and 
time-intermittent generation representative of possible fu- 
ture scenarios, we perform a centralized off-line stochas- 
tic optimization to select the generator-specific feedback 
parameters. These parameters need only be communi- 
cated to generators once per control period (60 minutes in 
our simulations). We show that inclusion of local power 



flows as feedback inputs is crucial and reduces frequency 
deviations by a factor of ten. We demonstrate our con- 
trol on a detailed transmission model of the Bonneville 
Power Administration (BPA). Our findings suggest that 
a smart automatic and distributed control, relying on ad- 
vanced off-line and system-wide computations commu- 
nicated to controlled generators infrequently, may be a 
viable control and communication architecture solution. 
This architecture is suitable for a future situation when 
generation-load imbalances are expected to grow because 
of increased penetration of time-intermittent generation. 

1 Problem Setup and Brief State- 
ment of Results 

In today's power systems, the system operator performs 
an Optimal Power Flow (OPF) dispatch periodically with 
typical time interval being 5, 15, or 60 minutes depending 
on the Balancing Area [1]. The OPF sets the power 
outputs of the committed generation to match power de- 
mand and minimize generation cost while respecting the 
capacity limits on lines, ramping constraints and limits 
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on generators and sometimes taking into account the N-l 
security constraints. In between two successive OPFs, 
the system is automatically controlled by a combination 
of two mechanisms. The faster of the two, acting on the 
scale of seconds, is primary frequency control-a fully 
distributed proportional feedback on locally-measured 
frequency deviations that may also include a deadband. 
The slower mechanism, acting on the scale of minutes, 
is automatic generation control (AGC), also called 
secondary control-a centralized feedback on the integral 
of a weighted sum of a centrally measured frequency and 
tie line flows to neighboring balancing areas.[2| 

These combined controls correct deviations in the 
generation-load balance driven by fluctuations in loads, 
renewables and other disturbances in the system. How- 
ever, these mechanisms do not explicitly incorporate 
line-flow limits, generators ramping limits, or time- 
integral constraints like those on run-of-river hydro 
generation or energy storage. For systems with relatively 
low levels of fluctuations, these limits are not frequently 
violated and it is not necessary incorporate them directly. 
However, higher levels of time-intermittent generation 
will create larger fluctuations and ramping events and 
the associated constraint violations will become more 
common. Standard primary and secondary controls are 
limited in their ability to balance these fluctuations, and 
better control design is needed to manage these larger 
fluctuations. Because these fluctuations are intimately 
connected to frequency deviations, they are of special 
concern because they my result in system-wide instabili- 
ties and loss of synchrony [3 1. 

Other considerations for real-time power grid control 
systems are communication constraints and commu- 
nication security[2]. Mechanisms that rely on central 
aggregation of the entire grid state followed by a centrally 
computed response will be vulnerable to communication 
failures and attacks on the communication network, 
making the overall system less robust. On the other hand, 
with significant renewable penetration, it is difficult to 
control a system purely based on local feedback, since 
under some conditions, it may be necessary to control 
distant generators in a correlated manner. 

In this preliminary work, we explore a hybrid approach 




Figure 1: Our model of the BPA Transmission Network 



that combines the speed and security of fully distributed 
control with the extensive system visibility provided by 
centralized control. Our method performs a centralized 
lookahead dispatch that also computes optimal local 
feedback parameters for all controllable generation, thus 
enabling the system to respond to fluctuations based only 
on local observables. We expand our definition of local 
observables to include not just frequency but also real 
power flows to neighboring nodes. We use an ensemble 
of forecasts that capture various possible scenarios for the 
wind generation and loads over the next intra-dispatch 
period (5 min/15 min/ 1 hour) to design an optimal 
time-varying dispatch for all the generators, as well as 
local feedback functions that enable the generators to 
respond to fluctuations based on the local observables. 

Our control design is split into 2 phases: 

a An off-line optimization phase where the distributed 
control gains are optimized jointly for the whole net- 
work in a central computer using extensive simu- 
lation of possible future wind generation and fore- 
cast scenarios. These gains are then communicated 
to each flexible resource (controllable device) in the 
transmission network. This off-line optimization 
would need to be re-run every time the statistics of 
possible future scenarios change significantly. In 
general, we expect this optimization to be run every 
time the generation re-dispatch changes. 

b An online response phase where each device imple- 
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ments its purely local control in response to local ob- 
servables (local frequency, line flows etc.) on the 
pace of the standard primary controls. 

We test our algorithm on historical data from the Bon- 
neville Power Administration (BPA) system |4|, an ideal 
test system for our algorithm as it has significant amounts 
of both hydro and wind generation. We show that our al- 
gorithm performs well, even in cases of significant wind 
ramps. 

Our results (detailed in section [3j> lead to the following 
important observations: 

a Local control based on response to frequency devi- 
ations and local line flows at each generation can 
keep frequency deviations down to the about 10 mHz 
while maintaining all the security and capacity con- 
straints. 

b Proportional control on frequency deviations and 
feedback on line flows is sufficient. Adding a 
frequency-deviation integral response is unneces- 
sary, which is advantageous because a distributed 
implementation of an integral term may cause insta- 
bilities due to errors in local frequency measurement, 
and also because it limits communication require- 
ments. 

c Joint optimization of feedback parameters for fre- 
quency deviation and line flows is necessary. Inde- 
pendent optimization or removal of either term leads 
to poor control performance. 

d Optimization over a finite but representative set of 
future scenarios enables the generalization of the 
control to new unseen scenarios. 

The rest of the paper is organized as follows: Section[2] 
describes the mathematical setting of the underlying con- 
trol/optimization problem; we describe and discuss results 
of our numerical BPA experiments in Section|3j and Sec- 
tion|4]presents conclusions and explains our path forward. 

2 Mathematical Formulation 

2.1 Preliminaries 

The power system is described by an undirected graph 
<$ = (S',y) with edges S and n vertices Y ' . The grid 



is composed of loads (1), conventional generators (g) and 
renewable generators (r). The flexible resources in our 
grid are the online conventional generators go. We denote 
by p the \Y\ x 1 vector of net active power injections at 
each node in the network and by p go ,p r ,p' the net active 
power injections due to online conventional generators, 
renewable generators and loads: p = p go + p r + p'. Each 
of these vectors is of size \Y\ x 1 with the convention that 
pf° = or p' = if there is no conventional or renewable 
generator at node i. Note here that we make the assump- 
tion that there is only one generator or load at a given 
node. If there are multiple, we replace them by an equiv- 
alent single generator or load. For a vector v with indices 
in Y , v, denotes a particular component for i € Y and \>s 
denotes the sub-vector {v,- : i € S} for a subset S C Y. 

2.2 Scenarios 

The control algorithm proceeds by analyzing an ensem- 
ble possible future scenarios and designs control strate- 
gies that optimize the system cost (defined below) across 
all scenarios. We define a scenario % to be a collection of 
the following quantities: 

a Renewable generation over the time horizon of inter- 
est: p r (f). 

b Load profile over the time horizon of interest: p|)(?). 

c A unit commitment (configuration of generators 
which are online, i.e. available for re-dispatch) go. 

To define the control problem, we require a collection of 
scenarios E and estimates for the probability of each sce- 
nario, i.e. E = ta,P r °b(X;)}- We note that for a given 
collection E, p r (f) and Pg(f) (items a and b from above) 
will vary across the ensemble of scenarios, however, we 
take go (the unit commitment from c) fixed because we 
are designing the time-dependent dispatch and local feed- 
back parameter for that particular go. In this work, we 
assume that the collection E is finite. Typically, E will be 
built up from load and wind forecasts from different fore- 
casting methodologies weighted by confidences in each 
of these forecasts. E could also include samples from 
a stochastic forecasting model based on climate models, 
historical data, meteorological sensors etc. 
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2.3 Control Formulation 

We ignore electro-mechanical dynamical transients and 
work with a discrete-time quasi-static approximation of 
the system dynamics with fixed time step 8 and integer 
time indices t = 0, 1, . . . , T: at each time step the power 
flows over lines are re-computed for configuration of 
consumption/generation at nodes evolving in discrete 
time. In general, the feedback can depend on any of 
the system variables, but we limit ourselves to local 
observables so that the control can be implemented 
in a completely distributed fashion at each generator 
after the dispatch and feedback parameters have been 
communicated. 

For each generator g <E go, we compute a time-varying 
dispatch p°(f) : < t < T, proportional frequency re- 
sponse coefficient a g , integral frequency response co- 
efficient a' g , and a response coefficient to local flows 
{(Xg^i ■ i € Neb(g)}. Further, we denote by 0)(t) the fre- 
quency deviation from the nominal frequency (50/60 Hz) 
at time t and by Cl(t) the integral of the frequency devia- 
tion, which in discrete-time is approximated by £l(t) = 
£' T=0 y T ~'(o(T) where < y < 1 is a discount factor. 
In other words, the integral frequency term is simply a 
weighted sum of frequency deviations in the past, where 
frequency deviations that are further in the past receive a 
geometrically smaller weight. With the time varying dis- 
patch and feedback parameters determined, the output of 
the generators is given by: 

Although our algorithm can incorporate nonlinear feed- 
back, we choose feedback which is linear in the local 
observables for this initial work. In addition to gener- 
ators, the real power consumption of loads responds to 
frequency changes, and we assume a simple linear load- 
frequency response given by 

P ! (0=Po(0 + P'fl>(0, 

where the j3' are known from measurement where pj, is 
the load at the nominal frequency (60 Hz). Combining the 
load and generator frequency response and the generators' 



time varying dispatch, the system's equilibrium frequency 
is computed by enforcing power balance in the system: 

LpT(0+p!(0+p'(0 = o =► 

, L-P!),(0+Pf(0+Pf (0 
LA- 

To compute power flow from the injections p(t) = p'(?) + 
p r (f ) +p g0 (f ), we use a modified version of the DC Power 
flow equations based on a linearization of the AC Power 
flow equations around the nominal dispatch at the begin- 
ning of the control period p(0). The linearization gives 
us dynamic impedances sf_^ ; that substitute for the line 
reactances in the DC power flow equations: 

p,(»=p»+ E 'M^, 

yeNeb(i) S i^j 

_ . frt __o | ft(0-Q>(0 
Vi^jl 1 ) -P i->7+ d 

Such a linearization is reasonable assuming that the flow 
patterns do not change too much during the course of the 
control period. 

In addition to several other constraints discussed below, 
we will also imposed a constraint on the total energy ex- 
tracted from generators in the control period. Such con- 
straints can represent the water discharge constraints on 
run-of-river hydro systems or state-of-charge constraints 
on energy storage devices. Therefore, we must also in- 
clude the total energy extracted from each generator into 
the system state: 

p'(0 = tp 80 ^)- 

T=0 



The overall system state consists of x(f) = 
[n(f);p 8 °(f);p 7 (f)] (in Matlab notation), and the 
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system evolution can be summarized by: 

1 ; LA 
Q(f + i) = a>(f) + yfl(0 

P g %(0 = P^W + «>W + «^(0 

r'eNeb(g) 

p I (0 = Po(0+JB i ®(0 

v 8/(0 n {t ,_ m-m 

yeNeb(i) *i->/ s i^y 

2.4 Cost Functions 

We consider a stochastic setting with many possible fea- 
tures, and it is unclear whether it is feasible to satisfy all 
constraints across all scenarios in E. Therefore, we use 
a penalty function to enforce our constraints in a smooth 
manner. The penalty function has a magnitude of zero in 
a dead-band around the most feasible region and grows 
cubically with the magnitude of constraint violation: 

Pen(a,/,M) = 

Here, a is the value of the constrained quantity and I 
and u are the lower and upper bounds on a, respectively. 
We also adopt the convention that when a,l,u can be 
vectors (of the same size) and the penalty in this case 
is applied element-wise and added up. The penalty 
function is designed so that the resulting cost function is 
smooth (twice differentiable). However, if a is violates 
the upper bound by 10%, a penalty of approximately 10 7 
is incurred-a high enough penalty so that if a feasible 
solution exists across all scenarios, it will be found. 

The cost function Cost(x(f),x(f + l),f) is computed at 
each time step in the control period, but it requires state 
information from both t and t + l so it can incorporate 
generator ramping limits. The cost includes seven terms 
that penalize both economic cost of supplying generation 
and deviations of the system state outside of normal oper- 
ational bounds. The individual terms are: 



1 Generation costs 

GenCost(p go (f)) = £ c g i(p| ) 2 +c g2 p| + ^ 3 . 

i'Ggo 

2 Generation limit penalties 

Pen(pS°( f ),p^,p^). 

3 Ramping limit penalties 

Pen(pr - 8 °" + '»-- 8 -"> .pf). 

4 Power flow thermal limit penalties 

E Pen(p<->X0.-P*->j.P)- 

5 Frequency deviation penalties 

Pen(fi)(f), -0.01,0.01) 

6 Integral frequency deviation penalties 

Pen(Q(f), -0.01,0.01). 

7 An integral deviation penalty on generation: 

Pen(p 7 (r),0.95£g°, 1.05£S°) 

Cost 1 simply represents the financial cost of energy from 
different generators. Costs 2-4 are normal power system 
constraints converted to costs using the penalty function 
defined above. Cost 5 is an additional penalty designed 
to constrain the system frequency to within a 10 mHz 
band, and Cost 6 is designed to constrain the deviation 
of the integral of the frequency deviation so that the fre- 
quency is not allowed to be low or high for extended pe- 
riods of time. Finally, Cost 7 is designed to keep the total 
energy delivered by each controllable generator over the 
control period within a ±5% band around a p 7 (r) mim- 
icking constraint the would occur in either a run-of-river 
hydro system or an energy storage device. 



'l0 7 ((a-M)/(0.1*( M +l))) 3 iffl> M 

< 10 7 ((/ -a)/(0.1 * (/ + 1))) 3 ifa</ 

otherwise 

V 
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2.5 Ensemble Optimal Control 



The evolution equations listed in ([T| are functions of a 
given scenario %, therefore, we can think of the state as a 
function of the scenario % and the control parameters a — 
{a p ,a',a F ,pf(t) : < t < T}: x{a,%,t). The overall 
optimization problem can then be written 



'T-\ 



min^Prob(j) £ Cost(x(a,x,t),x( a >X,t + l),t) 



,1=0 



Subject to ([TJ. 



(2) 



We optimize this objective using a standard numerical 
optimization algorithm (LBFGS 13). The gradients of 
the objective function can be computed efficiently using a 
forward propagation algorithm that uses the chain rule to 
propagate gradients in time. This computation can be eas- 
ily vectorized over all the scenarios, leading to significant 
speedup if run on a cluster or on GPUs. 




Time (minutes) 



(a) Aggregate Wind Generation 



V 



3 Numerical Results 



3.1 Description of Test System 

We test our algorithm using publicly available historical 
data for hydro and thermal generation, wind generation, 
and load from the Bonneville Power Administration 
(BPA) website||4). We use a model of the BPA trans- 
mission system (shown in Fig. [TJ that has 2209 buses 
and 2866 transmission lines. By identifying major 
hydroelectric stations on the transmission system and 
overlaying this onto a publicly available BPA wind site 
map ||6), we located the existing wind farms on the BPA 
transmission system (as of January 2010). We located 
the meteorological stations where BPA collects wind 
data Q in a similar manner. Using the same overlay, 
we used a simple incompressible air-flow model to infer 
hub height wind speeds at the wind farms. The resulting 
wind speeds were passed through the power curve of a 
standard 1.5-MW GE Wind Turbine which was scaled to 
the wind farm nameplate capacity to estimate the power 
output p r (?) (in MW) at each wind farm as a function 
of time. When we aggregate our wind farm-specific 
estimates of wind generation, we typically over estimate 
the BPA aggregate data by 20%, which may caused by 
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10 20 30 40 50 60 

Time (minutes) 

(b) Instantaneous and Integral Fre- 
quency Deviations 

Figure 2: Comparison of control schemes, a) Aggregate 
wind generation from a period with significant ramping 
events, b) Worst-case frequency deviations over the con- 
trol period for 1 8 validation scenarios not used in the con- 
trol design. 



6 



several factors including: spilling of wind by BPA, under 
performance of wind farms relative to single-turbine 
estimates, or shortcomings in our model of interpolating 
wind speeds. BPA also provides aggregate load data|4| 
that we divide among the nodes in the network according 
to population densities. BPA also makes publicly avail- 
able aggregate interchange flows [4|, which we apportion 
to different tie lines in a similar manner. 

To test our control algorithm on difficult conditions, 
we select a control period of one hour from 10:35 AM 
to 11:35 AM on February 12, 2010, when the wind 
generation was ramping significantly (shown in Fig. [2a| . 
We then create 26 scenarios (site-specific wind profiles) 
for this period by adding random time-varying Gaussian 
noise to the wind speeds at each meteorological station 
(from which we infer site-specific wind generation as 
outlined above). We set the magnitude of the noise so 
as to match, on average, the aggregate wind generation 
hour-ahead forecast errors reported by BPA [8|. All 
the time series data used in our study was available at a 
5-minute resolution. 

Unit commitment data is missing from our model, 
therefore, we assume that all hydro generators larger than 
300 MW are online and are all participating in frequency 
regulation. From inspection of the BPA historical genera- 
tion data [4], we infer that the thermal generation dispatch 
is fixed over time. In our model, we replicate this dispatch 
by dividing the total thermal generation among the online 
thermal generators (randomly chosen). 



3.2 Comparison 
Schemes 



of Various Control 



For difficult wind ramping conditions, we illustrate the 
value of feedback based on local flows by comparing four 
control schemes. We use P to designate proportional con- 
trol (to frequency deviations co(t)) and I designates inte- 
gral control (to integral frequency deviations £l(t)). The 
control schemes we consider using are: 

1 PI: Joint optimization of the time-varying dispatch 
Po°(0 and the local feedback parameters for co(t) 
and £l(f). 



2 Flow+PI Uncoordinated: Time-varying dispatch 
Po°(f) phis feedback on 0)(f),£2(f) and local flows 
Pg->; at each generator. The optimization in 1 is per- 
formed first followed by a second optimization over 
the flow feedback parameters. 

3 Flow+PI Coordinated: Same as 2, but the optimiza- 
tion is performed jointly. 

4 Flow+P: Same as 3, but without feedback on £l(t). 

The experimental protocol is as follows. We setup each 
of the four optimization problems according to Eqs. [2] 
with the scenarios described in Section [3~T1 and determine 
a single set of feedback parameters for each of the four 
feedback schemes. We use 8 of the 26 created scenarios 
as input to the optimization algorithm. The remaining 1 8 
unseen scenarios are reserved for validation of the control 
policy discovered by the optimization algorithm. We note 
that all four control strategies are able to achieve simi- 
lar generation costs while maintaining all the other con- 
straints (line thermal capacities, ramping limits, and inte- 
gral energy constraints), however, there are significant dif- 
ferences in the quality of the frequency regulation. Figure 
2b shows the worst-case frequency deviations over the 18 
validation scenarios. The frequency deviations are at an 
unacceptable level (.1-.2 Hz) when using just PI feedback 
(scheme 1). If the flow feedback is included but optimized 
separately (scheme 2), there is little improvement. How- 
ever, the if the PI and flow feedback are coordinated via 
joint optimization (scheme 3), the frequency deviations 
are reduced to an acceptable level. Interestingly, remov- 
ing the feedback on the integral of the frequency devia- 
tions (scheme 4) does not impact the frequency deviations 
significantly relative to scheme 3. 

3.3 Discussion of the Results 

The distributed frequency control method we have pre- 
sented benefits greatly from the incorporation of local 
power flows as demonstrated in Fig. [2b] There are several 
possible reasons for this improved performance. First, 
power flows make the local generation-load imbalances 
visible to the generators so that the closest generators re- 
spond, effectively screening the more distant generators 
from the need to respond. When compared to feedback 
based on frequency deviation, which is a global measure 
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of the imbalance, feedback on local power flows confines 
imbalances to shorter spatial scales with a corresponding 
decrease in the time scale of the response. An alternative 
explanation is that the optimization over the ensemble of 
possible futures in Eq. [2] is acting as a sort of machine 
learning that encodes correlations between the wind pre- 
diction errors and the resulting local power flows into the 
flow feedback parameters. When wind prediction error 
occurs, the change in power flows drives the feedback to 
nearly compensate for the error without a frequency devi- 
ation existing for any significant length of time. More nu- 
merical experiments are required to distinguish between 
these two (and other) possibilities. In both of the possibil- 
ities discussed above, variations in the local power flows 
appear to be acting as "pseudo-communication" channels 
between the renewable and controllable generators. Such 
a communication analogy may help explain why the in- 
dependent optimizations in scheme 2 does not yield sig- 
nificant improvement in control performance. The first 
optimization over frequency deviations may effectively 
washout the important local information in the power 
flows such that it is not available when optimizing over 
power flows. 

4 Conclusions and Future Work 

We introduced a control architecture based on off-line 
centralized optimization that can occur on a slow time 
scale coupled that sets the feedback parameters for fast 
distributed control of generation. The control scheme 
takes into account explicitly the variability in renewable 
generation using ensemble control. We showed that local 
feedback based on line flows and frequency deviations 
is sufficient to maintain all operational constraints and 
limit frequency deviations to an acceptable level even 
when the system is experiencing significant ramps in 
wind generation. Our method exploits the hour-scale 
predictability of wind energy while using the off-line 
optimization to re-adjust control policies over longer 
timescales where wind predictability suffers. Our hybrid 
approach has the potential enable even higher levels 
time-intermittent renewable generation than presented 
here, and it can do so without real-time computation or 
communication. 



These results are quite exciting and promising, how- 
ever, they are preliminary and much work needs to be 
done to ensure the viability of this scheme in practice. 

• Dynamical simulations are needed to check the dy- 
namical stability of a grid with flow feedback. If 
these simulations show that the scheme is unstable, 
we believe that this can be rectified by appropriate 
exciter control at the generators to damp the fast 
electro-mechanical transients. 

• The scenario approach can be extended to include the 
(N-l) security criterion, so that the optimized control 
strategy can deal with contingencies arising from the 
failure of a grid component. 

• It is possible that flow feedback acts as a pseudo- 
communication channel between generators in the 
absence of a dedicated communication channel. It 
would be interesting to investigate this from an infor- 
mation theoretic point of view and investigate how 
much of information can be encoded in the flows. 

• We have used the simplest possible algorithmic ap- 
proach by defining a smooth version of the optimiza- 
tion problem using penalty functions solving it us- 
ing a generic LBFGS algorithm f5). Second-order 
algorithms such as Stagewise NewtonJ3 or Differ- 
ential Dynamic Programming (DDP)[10| efficiently 
exploit the problem structure of deterministic opti- 
mal control problems. These can be leveraged in our 
ensemble control context by noting that when the 
feedback parameters a' ,a p \a F are fixed, we have 
a deterministic optimal control problem in Po°(f) 
for each scenario. We have also been working on 
a Gauss-Newton algorithm for optimizing the fixed 
feedbacka 7 , a p , a F efficiently. One can perform al- 
ternate minimization of pf?(t) and a 1 , a p , a F to get 
an efficient algorithm for optimizing both. Further, 
we note that when feedback does not include the in- 
tegral term £i(f), the ensemble control problem is a 
convex programming problem, and the global opti- 
mum can be found efficiently using specialized con- 
vex optimization techniques. 

• We plan to incorporate more accurate AC model- 
ing of power flows taking advantage of most recent 
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advances in analysis and algorithms related to opti- 
mizations of nonlinear power flows, e.g. lfTTl[T2ll . 

• The integral energy constraint we introduced can 
also model energy storage, and our algorithm can 
easily be extended to incorporate distributed control 
of energy storage. 
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